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Inspired by "quantum graphity" models for spacetime, a statistical model of graphs is proposed to 
explore possible realizations of emergent manifolds. Graphs with given numbers of vertices and edges 
are considered, governed by a very general Hamiltonian that merely favors graphs with near-constant 
valency and local rotational symmetry. The ratio of vertices to edges controls the dimensionality of 
the emergent manifold. The model is simulated numerically in the canonical ensemble for a given 
vertex to edge ratio, where it is found that the low energy states are almost triangulations of two 
dimensional manifolds. The resulting manifold shows topological "handles" and surface intersections 
in a higher embedding space as well as non-trivial fractal dimension. The transition is first order, 
underlying the difficulty of graph models in describing criticality that is independent of the details 
of the underlying graph. Another interesting phenomenon is that the entropy of the graphs are 
super-extensive, a fact known since Erdos, which results in a transition temperature of zero in the 
limit of infinite system size: infinite manifolds are always disordered. Aside from a finite universe or 
diverging coupling constraints as possible solutions to this problem, long-range interactions between 
vertex defects also resolve the problem and restore a non-zero transition temperature, in a manner 
similar to that in low-dimensional condensed-matter systems. 



I. INTRODUCTION 

Since the first systematic studies of random graph 
models by Erdos and Renyi [T], the relation between 
graph theory models and physics models, in particu- 
lar statistical physics models, has attracted much inter- 
est. Concepts and tools in graph theory have been ap- 
plied to problems in physics, computer science, and biol- 
ogy to produce remarkable results. For example, Feyn- 
man diagrams which are planar have special roles in the 
large N QCD model P]; in causal dynamical triangu- 
lation, 4-dimensional triangulated manifolds with fixed 
edge lengths, which can be viewed as a class of graphs, are 
used to construct spacetime on the Planck scale to regu- 
larize the quantum gravitational path integral [31 3] ; sta- 
tistical mechanical models of network growth can explain 
the connectivity of systems such as the internet [5 ; struc- 
tures of amorphous solids can be quantified using graph 
theory properties [6]; intracellular signalling networks 
can exhibit emergent behavior stored within biochemical 
reactions, including integration of signals across multiple 
time scales and self-sustaining feedback loops [7]; neural 
networks can collectively and robustly produce content- 
addressable memories from partial cues [8 , indicating 
capacity for generalization, familiarity recognition, and 
categorization. Added to these discoveries, a new collec- 
tion of graph models have been proposed as candidates 



for emergent spacetime, as described below. 

A manifold can be approximated by a triangulation, 
which in turn can be viewed as a graph filled with sim- 
plices. From this observation, one can consider how a 
graph may give rise to a manifold, i.e. from a family 
of graphs, following some constraints and obeying some 
set of rules for dynamical processes, is it possible that 
a manifold-like structure can emerge? To be more pre- 
cise, consider the possibility that a graph G gives rise 
to a smooth manifold M. A vertex in G corresponds 
to a point in AI; when a pair of vertices in G are con- 
nected by an edge, the corresponding pair of points in 
M have a certain distance e. When the length scale un- 
der consideration is much larger than e, G resembles the 
smooth manifold M . In such cases, one can say that the 
manifold M , including its dimensionality, topology, and 
metric, emerges from the graph G in the continuous limit. 

From this general idea, in references [IJ [10] , a graph 
model was constructed from a given graph Hamiltonian, 
where it was proposed that the low energy phase of the 
model may be interpreted as an emergent spacetime. In 
addition, it was found that when the edges of the graph 
possess a spin degree of freedom, the model could give rise 
to a U(l) gauge theory [TO]. In [TT], Konopka has ana- 
lytically and numerically studied the above graph model 
as a statistical model. A phase transition was found, 
where it was argued that the low temperature phase can 
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be related to spacetime only if the graph can interact 
with some matter degrees of freedom. In [T^ [T3], a re- 
lated model, which in addition to graphs corresponding 
to spacetimes, also incorporates a matter field which re- 
side on the vertices, is proposed to study the role of 
matter in the emergence of spacetime from graphs. In 
[14j . Conrady has constructed a Hamiltonian favoring 
low temperature, 2-dimensional manifolds through terms 
that explicitly favor 2-d triangulations, for example each 
vertex is favored to have 6 edges as in a triangular lattice, 
and tetrahedra are penalized. The model was simulated 
for small system sizes {N < 180 edges), which showed 
a heat capacity peak, and a transition temperature that 
decreased with system size. 

In this paper, we also investigate a statistical model 
of graphs, in that the objects under consideration are 
merely abstract graphs, without any information on the 
positions of the vertices, or the lengths of the edges. A 
graph can randomly transform into another graph ac- 
cording to a set of transformation rules. Graphs with 
given numbers of vertices and edges are considered, and 
they are governed by a Hamiltonian that favors graphs 
with a set of local symmetries. If these local symme- 
tries are preserved, the resulting graphs should be nearly 
triangulations of manifolds with a certain dimensional- 
ity, where the dimensionality is controlled by the ratio 
of vertices to edges. We are interested in whether any 
global structure of the graphs arises as a consequence. 

Because every edge in this model corresponds to a posi- 
tive length e, only real positive distances can arise, so this 
model can only be used to describe Riemannian mani- 
folds (i.e., with positive definite metric). The metric of 
a Riemannian manifold can be alternatively viewed as a 
distance function between any pair of points, which sat- 
isfies the triangle inequality. On a graph, there is also 
a natural notion of distance, namely the length of the 
shortest path between a pair of vertices. This distance is 
also positive-definite and satisfies the triangle inequality. 
Thus on any 

graph, there is a well-defined distance function, as well 
as a corresponding geodesic. Graph geodesies between 
two vertices are often highly degenerate however, unlike 
the case for manifolds. If a manifold is to emerge from a 
graph, one expects that in the continuous limit all degen- 
erate geodesies are close by, and the differences of their 
path lengths are only of order e. After establishing this 
distance function between vertices, mapping the graph 
to a Riemannian manifold is still a non-trivial problem. 




FIG. 1. Examples for the graph theory concepts. 



If we enforce that every edge is identical in that they 
have the same length when mapped to the Riemannian 
manifold, then only for certain graph configurations will 
a Riemannian manifold emerge from the graph. Other- 
wise the system will be frustrated and unable to meet the 
condition of constant edge length, without increasing the 
dimension above that of the manifold that would emerge 
from the graph. 

In this paper, after reviewing the relevant graph the- 
ory preliminaries, we introduce a graph Hamiltonian 
based only upon local symmetries. We evolve the graph 
under the Monte Carlo rules obeying statistical me- 
chanical equilibrium, and we investigate whether a low- 
temperature manifold state emerges. We investigate the 
sharpness of the phase transition using energy as an or- 
der parameter for different size systems, and discuss the 
likely first-order nature of the transition. We construct 
heat capacity curves as a function of temperature, and 
investigate the transition temperature as a function of 
system size, which points towards a zero-temperature 
phase transition in the bulk limit. The Haussdorf di- 
mensionality of the emergent manifold is investigated, 
and found to be an increasing function of system size, 
and approximately 3 for the largest system sizes we in- 
vestigated (2000 vertices). Correlation functions between 
defect-carrying vertices and edges are investigated to de- 
termine whether the effective potential between defects 
is attractive or repulsive. Lastly, we argue in analogy to 
condensed-matter systems that a non-zero phase tran- 
sition temperature requires long-range interactions, and 
show that a Coulombic-likc term between graph vertices 
yields an apparently finite-phase transition temperature, 
but with a highly ramified manifold. 



II. GRAPH THEORY PRELIMINARIES 

Before motivating for details of the model, we shall 
remind the reader about some graph theory concepts, 
which will be needed later in constructing the model. 

A graph G is composed of a set of vertices V{G), 
and a set of edges E{G), where every edge is a subset of 
V{G) with two elements. Note that by this definition, 
the two vertices in an edge set cannot be the same ver- 
tex, and two edges cannot connect the same two vertices. 
Such graphs are sometimes called "simple graphs" , as op- 
posed to "multigraphs" . Because we will only consider 
graphs of this definition, they will be simply referred to 
as "graphs". 

A vertex v is incident with an edge e ii v G e. We 
denote an edge e by its vertices, or ends, say u and v, as 
e = {u,v}, or simply e = uv. A vertex m is a neighbor 
of, or is adjacent to, a vertex v if uv is an edge. The 
valency or degree of a vertex is the number of edges 
incident to that vertex. 

A graph in which every vertex has the same valency is 
regular. It is fc-regular if every vertex has valency k. 

A graph in which every pair of vertices is connected 
by an edge is complete. It is denoted by Kn if it has n 
vertices. 

G' is a subgraph of a graph G is a graph if V^(G') C 
V{G) and E{G') C E{G), and this is denoted by G" C G. 

If L/ C V{G), the subgraph G' induced by U is the 
graph for which V{G') = U, and E{G') contains an edge 
xy if and only if x, y e U and xy G E{G). This is denoted 
by G' — G[U], and G' is called an induced subgraph of 
G. (For example, in Figure [l] the vertices k,o,p,s, and 
the five thick edges, compose an induced the subgraph; 
the vertices i,m,n,q^ and the four thick dotted edges, 
compose a subgraph, but not an induced subgraph.) In 
particular, in a graph G, the subgraph induced by the set 
of neighbors of a vertex v is called the neighborhood of 
V, and is denoted by Gn{v). 

A path is an alternating sequence of vertices and 
edges, beginning with a vertex and ending with a ver- 
tex, where each vertex is incident to both the edge that 
precedes it and the edge that follows it in the sequence, 
and where the vertices that precede and follow an edge 
are the end vertices of that edge. The length of a path 
is the number of edges in the path. (For example, in 
Figure [l] (a, ab, b, bf, /, fg,g) is a path with length 3, in 
which the edges are denoted by dotted lines, and is also 
one of several paths between a and g having the mini- 




FIG. 2. Some examples of neighborhood subgraphs. The ec- 
centricity is labeled for each vertex, and the difference of di- 
ameter and radius of these subgraphs, which is denoted by A, 
is labeled below each graph, (a) is the neighborhood subgraph 
of vertices in the triangular lattice graph (Figure [T]); (b), (c), 
(d) and (e) appear commonly in simulations, as parts of the 
defects; (f) is the neighborhood subgraph of vertices in the 
graph in Figure |3] 



mal distance.) The distance between two vertices is the 
length of shortest path between them. In a graph G, the 
distance between vertices u, v is denoted by da{u, v). 

A graph is connected if any two vertices are linked 
by a path. 

The eccentricity eciv) of a vertex u in a graph G is 
the maximum distance from v to any other vertex, i.e., 

eoiv) = max daiv^u), 

u£V(G) 

where dciv^u) is the distance between v and u in the 
graph G. 

The diameter diam(G) of a graph G is the maximum 
eccentricity over all vertices in a graph, and the radius 
rad(G) is the minimum, 

diam(G) — max eniv), rad(G) — min enM- 

vev{G} veV{G) 

When G is not connected, diam(G) and rad(G) are de- 
fined to be infinite. Some example of neighborhood sub- 
graphs are shown in Figure [2j For every vertex in Figure 
[T] the neighborhood subgraph is Figure ; for every 
vertex in Figure [3j the neighborhood subgraph is Fig- 
ure [2]jf). Figures [2j^b)-[2f^e) are examples of neighborhood 
subgraphs that appear commonly in the simulation. 

Given a lattice, the corresponding lattice graph is 
the graph whose vertices are the points in the lattice, 
and whose edges are the pairs of nearest points in the 
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lattice. (For example, the whole graph in Figure [T] is an 
equilateral triangular lattice graph.) 

III. THE MODEL 

To gain intuition on the form of constraints and Hamil- 
tonians that may induce manifolds, let us construct some 
graphs resembling some manifolds, starting with the ex- 
ample of a flat two-dimensional plane . Intuitively, any 
two dimensional lattice graph as defined above forms a 
"two dimensional" manifold, and a coordinate system of 
the manifold naturally inherit from the coordinates of 
the lattice graph. This is directly analogous to a Bra- 
vais lattice in crystallography. A priori there seems no 
decisive reason to choose any particular Bravais lattice 
as the preferred graph configuration, however we shall 
choose the equilateral triangular lattice graph (Figure 
[l]), using the following argument. On E^, for any point 
p and any distance S, let Bs{p) denote the geodesic ball 
centered at p with radius 6, and Bs{p) —p has the topol- 
ogy of a circle S^. For graphs, we can define the notion 
of "geodesic ball" similarly with that in Riemannian ge- 
ometry. Let Bn{v) be the set of the vertices that have 
distance from vertex v no greater than n, including v 
itself. For any two dimensional lattice, if we denote the 
corresponding graph by G, for sufficiently large n, the in- 
duced subgraph G[i?„(w) —v] also looks like topologi- 
cally. However, for n = 1, namely the neighborhood sub- 
graph Gn{v) = G[Bi{v) — v], this property is no longer 
true for all lattices. For example, on the square lattice, 
Gn{v) is composed of 4 disconnected vertices. Only for 
the equilateral triangular lattice, G]y{v) looks topologi- 
cally like S^. Thus in this sense, equilateral triangular 
lattice graph is the closest analog to among all the two 
dimensional lattice graphs, on all distance scales down 
to e. 

A graph can form a two dimensional lattice for the cor- 
rect ratio of edges to vertices. While a thermalized lattice 
in two dimensions is isotropic |15H17| , the connectivity of 
such a lattice is still well-defined at low temperature. We 
thus choose to add defects in the form of extra edges or 
bonds, which will evolve under some Hamiltonian. This 
allows bonded vertices to be permuted, so that the low 
temperature phase is still a "quasi-fluid" that retains a 
symmetry corresponding to randomized graph connectiv- 
ities. The extra edges induce defects in the lattice, which 
may be mobile. The exact shape of the defects and the 
reason why the defects are unstable or meta-stable de- 



pend sensitively on the Hamiltonian. We shall construct 
a candidate Hamiltonian, and test the stability of the de- 
fects by numerical simulation. This construction general- 
izes to ]R" straightforwardly: We can see that the defect- 
free lattice is the rt-dimensional lattice as a arising from 
a regular tiling of n-diniensional tetrahedra. The defect 
is a (n — l)-dimensional "foam" that divides the space 
into many patches of lattices with random orientations. 

We seek the simplest Hamiltonian that can give rise to 
manifold-like triangulation graphs as classical solutions, 
which contain defects that facilitate graph permutation 
symmetry. We assume that the action is local, in the 
sense that it should be a sum over the vertices and/or 
edges, such that each term involves a finite number of 
vertices and/or edges within some cutoff distance. This 
condition is imposed because almost all physics models 
for which the Hamiltonian or Lagrangian is an integral 
of the corresponding density are local in the same sense. 

A defect manifests itself as a local structure contain- 
ing vertices with anomalous valency. One obvious local 
property of manifold-like graphs is that all vertices not 
in any defects would have the same valency. Moreover it 
is likely that vertices in the defects have just one more 
or one less neighbors. These properties can be enforced 
by a Hamiltonian quadratic in the valency: 

vev{G) 

where n„ is the valency of vertex v, and Ci is a positive 
constant (which will be taken to be infinite as described 
below). The average valency of the vertices is given by 



where Ne is the total number of edges and Ny the total 
number of vertices. Note that for example a = 6 is com- 
patible with a regular equilateral triangular lattice, which 
in turn implies that the emergent manifold is two dimen- 
sional, while a = 12 is compatible with the face-centered 
cubic lattice, which implies a three dimensional emergent 
manifold. Thus without changing the form of the Hamil- 
tonian, we should be able to find manifolds with different 
dimensionalities by adopting different a priori values of 
a. In the simulations described below, we choose a to 
be a non-integer, so that there exists an "excess" num- 
ber of edges, which contribute to the presence of defects. 
Because the total number of vertices and edges are fixed, 
the term in ([I]) is minimized when every vertex has va- 
lency either [aj or \a \ . In our simulations, ci is taken to 
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FIG. 3. A 6-regular graph which is not similar to any man- 
ifold. This graph can be viewed as an infinite rooted "tree" 
graph, in which each node has three children (except the root 
node has four children), and every node of the "tree" is actu- 
ally a tetrahedron. 

be infinite and so is no longer an adjustable parameter of 
the model, and the corresponding term in ([T]) is enforced 
to be minimal. 

To obtain manifold-like solutions consisting of patches 
of close-packed lattices interspersed with defects, it is not 
sufficient to only impose the condition that each vertex 
has approximately the same number of neighbors. Many 
regular graphs do not look like any manifold at all (see for 
example Figure |3| . Additional terms in the Hamiltonian 
are thus required for manifold-like solutions. 

One candidate for such a term consists of particular 
subgraphs that can be embedded into the graph. From 
this viewpoint, n„ is the number of K2 subgraphs (two 
vertices connected by an edge) that go through the vertex 
V. It is likely however that choosing more terms of this 
type will affect the dimensionality of the resulting space- 
time. For example, if we incorporate terms that favor 
more K3 subgraphs (triangles) and less K4 subgraphs 
(tetrahedra) , then it can be expected that these terms 
would favor two dimensional manifolds |14j . As we hope 
to find a model which does not select the dimensionality 
at the level of the Hamiltonian, we will not use any other 
term of this type besides Hi . 

Another property of manifold-like graphs is that 
around most vertices, the graph has a local discrete ro- 
tational symmetry that reflects the local isotropy of the 
emergent manifold. This can be restated as for each ver- 
tex t), the subgraphs G[B„(u) — v] for most v should have 
a discrete rotational symmetry. To reduce the number of 
possible Hamiltonian terms, we only impose this condi- 
tion on G[Bi{v) — v], which is also Gn{v). We introduce 
the term 

H2 = C2 A(i.)' (3) 

vev{G) 



where C2 is a positive constant, and 

A{v) = diam(Gjv(w)) - rad(GAr(t;)), (4) 

in which diam(G ^ (v)) is the diameter of the subgraph 
Gn(v)^ and rad(Gjv(f)) is the radius of the subgraph 
G]sf{v). By the definitions of diameter and radius of 
graphs, if the subgraph Gpf{v) is not connected, they are 
both infinite. Here, we additionally define that their dif- 
ference diam(GAr('(;)) — Tad{G n (v)) is also infinite when 
Gn(v) is not connected. The term H2 then enforces that 
all neighborhood subgraphs are always connected. When 
Gn{v) is connected, the difference between its diameter 
and its radius is a measure of its asymmetry. Figure [2] 
shows several examples of neighborhood subgraphs. The 
eccentricity of every vertex in the subgraphs is labeled, 
along with the value of A{v) for each subgraph. For Fig- 
ures [2ja) and (b), the GAr(i')'s have a rotation symmetry 
of Dg and D7, respectively, while Figures [2jc)-(e) are not 
rotationally symmetric. 

In two dimensions, a graph forms a triangulation of a 
surface if and only if all the neighborhood subgraphs are 
cycles [H] . When the degrees of the subgraphs are either 
6 or 7, which is imposed by the Hi term, one can see from 
the examples in Figure [2] that the H2 term indeed favors 
cyclic neighborhood subgraphs, with only one exception 
shown in Figure[2je). We thus expect that, in this model, 
a graph with low energy is almost a triangulation of a 
surface. 

Thus we propose the following model: Consider a sim- 
ple graph with Ny vertices and Ne edges. All the ver- 
tices are labeled, so isomorphic configurations with differ- 
ent labeling are considered to be different configurations. 
The Hamiltonian is composed of two terms, as motivated 
previously: 

H = Hi+H2. (5) 

Because the Hamiltonian is prohibitive to analytical 
solution, we implement a numerical simulation, as de- 
scribed in the next section, to study the equilibrium 
states of this model in the canonical ensemble, i.e. at 
a given temperature. In particular, we will be interested 
in the structures of the states with low energies, and the 
nature of the phase transition, if one exists, to these low 
energy states. 
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IV. NUMERICAL SIMULATION 

We sample equilibrium states in the model using a 
Monte Carlo simulation [19]. The parameter a defined 
in ([2]) as giving the mean number of edges per vertex is 
taken to be slightly larger than 6, which we expect will 
induce two-dimensional structures dictated by triangu- 
lations as described above. There is no fixed boundary 
on the graphs. The size of the graphs is specified by the 
number of vertices Ny, and the number of edges Ne- For 
convenience, in the following we use Ny and the number 
of extra edges X = Ne — ^Ny , to specify the size of the 
graphs. Given the graph size, the initial configuration is 
taken to be a randomly generated, connected graph. 

The graph is evolved in the canonical ensemble with 
temperature 1/(3. In each Monte Carlo step, one end of 
an edge can jump from one vertex to another. We ran- 
domly pick an edge, and randomly label its ends by u and 
V. To find the new location of the edge uv, we perform 
a random walk starting from v as the origin, which does 
not pass through the edge uv (this condition guarantees 
that a connected graph remains connected after such a 
move). The number of steps £ of the walk is a random 
positive integer chosen from the probability distribution 
P{£) — J^^^ — 7*^, where 7 is a parameter between and 
1 (we take 7 = 0.5 below). Denote the ending vertex of 
the random walk as v' . The edge is then moved from 
uv to uv' . If the new graph is still simple, its energy 
is compared with that of the old graph, and this move 
is accepted or rejected according to the Metropolis algo- 
rithm [19J. Each "sweep" through the system contains 
Ne Monte Carlo steps, so on average each edge has one 
chance to jump in one sweep. Such a method is ergodic; 
moreover with this jumping scheme, the energy of only a 
few vertices is affected after each Monte Carlo step, and 
the energy of only these vertices needs to be updated. 

Simulations are performed with ci = 00, C2 = 1.0, 7 — 
0.5, and various values of Ny, X and /3. Before showing 
the thermodynamics results from the simulations, let us 
first describe the method that we used to render a graph 
from the simulations, in order to interpret its evolution. 



A. Rendering graphs 

To render a graph such that its structure can be best 
visualized, we need to devise an appropriate drawing 
scheme. A drawing of a graph maps vertices to points in 
M" with line segments connecting adjacent points. The 



following method is used to generate drawings in M^. For 
any drawing of a graph G, we seek to minimize the func- 
tion 



E 

e£E{G) 
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i^j not adjacent 



(6) 



where le is the length of the drawing of edge e, Uj is 
the distance of the drawing between vertices and 
fli = 1.0,02 — l-Oi^s = 5.0. The first term gives a pre- 
ferred length for every edge, and the second term gives 
a repelling force to every non-adjacent pair of vertices. 
The function i?draw is chosen this way in order to make 
every edge have approximately the same length in the 
drawing, and as well, to make the drawing as expanded 
as possible. In practice, even for moderate-sized graphs, 
Hdvam has numerous local minima and is difficult to min- 
imize. We thus use another Monte Carlo calculation to 
search for its near-optimal values. Initially, all the ver- 
tices are located at the origin of M^. In each Monte Carlo 
step, a randomly-chosen vertex is randomly moved to an- 
other position within the ball of radius 5 = 2.5, centered 
at the original position, and the new position has uni- 
form probability distribution within the ball. After the 
Monte Carlo calculation, because the low temperature 
configurations in the model are conjectured to be similar 
to triangulations of surfaces, we also search for all the 
subgraphs (triangles) in the graph, and render (flat, 
solid) triangles to fill the interior of the K^^s. 

Figure |4] shows some snapshots taken from the sim- 
ulations. Figure |4|^a)-|4][b) are for the system of size 
Ny 200 and X = 20. Figure [Ij^a) shows the initial 
configuration, |4][b) shows a typical configuration at high 
temperature (/3 = 1.0), and|4]jc) shows a typical config- 
uration at low temperature (/3 = 2.0). Figure |4]jd) is for 
the system of size Ny — 1000 and X = 100, and it is a 
typical configuration at low temperature (/3 = 2.0). 

In the sample drawings in Figure |4j different colors are 
used to denoted different types of vertices. The color- 
code is as follows: 



zero contribution to H2 
nonzero contribution to H2 



degree^ 6 
black 
red 



degree^ 7 
green 
blue 



Also, yellow lines are drawn at places where two trian- 
gles intersect, i.e., this identifies where the triangulated 
surface intersects with itself. 



(a) 



(b) 





(c) 



(d) 



FIG. 4. (Color online) Some snapshots from the simulations, drawn in three dimensions. Panels (a)-(c) are for the system 
with number of vertices Ny ~ 200 and number of extra edges X — 20, where (a) is the initial configuration, (b) is a typical 
configuration at high temperature (/3 — 1.0), and (c) is a typical configuration at low temperature (/3 = 2.0). Compared with 
the sphere, the drawing (c) has three more handles, and the surface intersects with itself in three places, so it has a non-trivial, 
non-orientable topology. Panel (d) is for the system of size A*'^ = 1000 and X = 100, and shows a typical configuration at low 
temperature (/3 = 2.0). In these drawings, if a vertex has valency 6, it is black if its A value is zero, and is red if its A value 
is nonzero; if a vertex has valency 7, it is green if its A value is zero, and is blue if its A value is nonzero (see text). As well, 
yellow lines are drawn at places where two triangles intersect, and the manifold thus passes through itself. 



B. Topology of the manifold in the presence of 
defects 

For the low temperature graphs, several examples of 
common local defects are shown in Figure [Sj They arc 
called local in the sense that in the vicinities of these 
defects, the graph is similar with some triangulation of 
surfaces with trivial topology. Among these examples, 
(a)-(c) do not increase the total energy, and around such 



defects the ratio between the number of edges and ver- 
tices is larger than 3. In other words, these defects can 
"absorb" the extra edges without energy cost. Also note 
that (a) and (b) do not change the long range order of 
the lattice orientation, while (c) does alter the long range 
order. Taken together, these defects induce configura- 
tional degeneracies in all the energy levels, including the 
ground state energy level, and at the same time induce 
graph permutation symmetry by randomly breaking the 
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FIG. 5. Example of some common defects. The figures in the first row are the schematic drawing the defect, in which a vertex 
is marked with a square if its valency is 7, a vertex is marked with an open circle if it contributes positive energy to and 
otherwise a vertex is marked with a filled circle. The figures in the second row are the corresponding drawings of the defects 
using the method described in subsection IV A Compared with the equilateral triangular lattice, examples (a), (b), (c), (d) 
and (e) have 2, 3, 2, 0, extra edges, respectively. 



lattice's long range order, at least in the rendering scheme 
of the manifold described above. Defects (d) and (e) in- 
crease the total energy, and alter the lattice orientation 
more drastically. 

As discussed above, low-temperature graphs in the 
model are similar to two dimensional triangulated sur- 
faces. However, they contain local defects, and there are 
overall topological features of the surfaces that emerge 
from the graphs. For example, in the drawing Figure 
^c) , one can see that the emergent surface contains sev- 
eral handles, and the surface intersects itself in several 
places. In the drawing Figure Qd), the topology of 
the emergent surface is too intricate to easily identify. 
The Hamiltonian does not constrain the topology in any 
way, so in general, emergent surfaces of low-temperature 
graphs in the model have complicated topologies. The 
emergent surfaces have potentially many handles, and 
are in general non-orientable, in that there is no separa- 
tion between interior and exterior sides of the surface. In 
our simulations, we also observe that the topology of the 
graphs' emergent surfaces can dynamically change, even 
at a low energy. 

We note however that the choice of Ny and Ne can 
constrain the topology. At low temperatures, the graphs 
are nearly triangulations, albeit with potentially com- 
plicated topologies. If a graph is strictly a triangu- 
lation, and we denote the number of triangles as Np, 
then the Euler characteristic x of the surface is given by 
^^ = Ny - Ne + Np- For a triangulation, SNp = 2Ne; 



and we previously defined Ne — 3iVy + X. Putting 
these three equations together, we find x = —X/3. As we 
showed above, defects on the graphs can absorb edges, so 
the relation for the nearly-triangulated graphs becomes 
an inequality x ^ —X/3. In addition, for any surface, 
X < 2, with X = 2 corresponding to the topology of a 
sphere. Thus the Euler characteristic x of the emergent 
surface can take any integer value between —X/S and 2. 
The X values used in our simulations are not very small, 
so this constraint still allows for many possible different 
topologies for the emergent surface. 

C. Phase transition 

In this sub-section we study the transition between 
the low-/high- temperature phases. For system sizes 
Nv = 100, 200, 300, 500, 1000, 1500, 2000, and number of 
excess edges X = O.lA^y, the expectation value of en- 
ergy (E), and the heat capacity C = P'^ {{E'^) - (E)'^) 
are computed for various inverse temperatures /3, where 
the angle bracket here means averaging over all the sam- 
ples in a simulation. 

The results are shown in Figure [6] and Figure [7j For 
the three largest systems with Ny = 1000, Ny = 1500 
and Ny — 2000, we also employ the weighted histogram 
analysis method (WHAM)|20) [5T] to improve the sam- 
pling quality. The inverse transition temperature /3c is 
defined as the inverse temperature where the heat capac- 
ity is maximal. It can be seen that /3c increases as Ny 
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FIG. 6. (Color online) The average energy density (E) /Nv 
as a function of inverse temperature /3 for for several TVv's 
indicated in the legend. 



FIG. 8. Log-log plot of the inverse transition temperature jSc 
in the model as a function of system size Nv, and the best 
fit line. The straight line fit indicates that as Nv — > oo, the 
transition temperature Tc — >■ 0. 
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FIG. 7. (Color online) The rescaled heat capacity C /Nv = 
P'^ {{E"^) - {Ef) /N^ as a function of inverse temperature /3 
for several A^v's indicated in the legend. 



increases. Near the transition temperature /3c, d{E) /df3 
also increases as Ny increases, and thus the widths of the 
heat capacity peaks decrease as Ny increases, indicating 
the transition becomes more cooperative. Figure |8]shovifS 
a log-log plot of the inverse transition temperature as a 
function of Ny ■ The linear relation in the plot indicates 
that as Ny goes to infinity, the transition temperature 
would go to zero. In addition. Figure [9] shows the prob- 
ability density distribution of E/Ny, for the systems of 
size Ny — 1000, 1500 and 2000, at each system's tran- 
sition temperature. As Ny increases, the energy distri- 
bution of the two phases become more bimodal, which, 
along with the trend in cooperativity with system size in 
Figure |6] indicates that the transition is first order in the 
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FIG. 9. (Color online) The probability density of the intensive 
energy E/Nv for the systems of size A'^v' = 1000, 1500 and 
2000, at each system's transition temperature. The error of 
p{E/Nv) for E/Nv > 0.5 is small (Ap < 0.1), the error for 
0.01 < E/Nv < 0.5 is Ap < 0.6, and the error for the smallest 
values of energy E/Ny < 0.01 is Ap < 2.5. 



bulk limit, with a corresponding nucleation barrier. 

In appendix|B] we give the acceptance ratio in our sim- 
ulations as a function of inverse temperature. Although 
the acceptance ratio substantially decreases in the low 
energy phase, the system is still able to undergo dynam- 
ics because some local defects cost little or no energy. 

A transition temperature of zero for infinitely large 
graphs is actually not very surprising on entropic 
grounds. Consider a first order phase transition of an 
extended physics model. Denote the size of the system 
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FIG. 10. The entropy density difference across the transition 
AS/Nv as a function of Ny- The best fit line using a loga- 
rithmic function is also shown. The inset shows AS/Ny as a 
function of Ny for a model including a Coulomb potential be- 
tween valency-7 vertices (see section |v|. Including long-range 
interactions can remove super-extensivity of the entropy. 



by N, and denote the number of states in the high and 
low temperature phases by JIh and f^L, respectively. Be- 
cause the energy difference between these two phases is 
proportional to N, the phase transition temperature Tc 
is given approximately by riue^"^^^'' = flh, where k is a 
positive number. As N increases, for a "typical" physics 
system with short-ranged interactions, the ratio between 
r^H and r^L increases as e"'^ , where 7 is a positive num- 
ber. This behavior results in a finite, non-zero transition 
temperature in the infinite size limit. On the other hand, 
the number of inequivalent graphs with Ny vertices is 
typically Ny , (for example, see [Ullll]) where 7' is a 
positive number that depends on the constraints of the 
allowed graphs. In our case, the allowed graphs should 
have every vertex valency equal to six or seven, and ev- 
ery vertex neighborhood should be connected. While we 
do not have an algorithm to count the exact number of 
allowed graphs, it is reasonable to assume for our sys- 
tem that the ratio between fiu and J^l has the typical 
asymptotic behavior of graphs, which explains a transi- 
tion temperature of zero, i.e., the transition temperature 
Tc is given by N^'^'' e"''^''/^'' w 1. 

To validate the above argument, we can calculate the 
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FIG. 11. (Color online) Distribution of edge valencies as a 
function of inverse temperature /3c, for the system of size 
Nv ~ 1000. There are no edges in the simulation with edge 
valency less than one or larger than five. 

entropy difference across the transition as given by 

where C is the heat capacity, and /3i and P2 are typical 
inverse temperatures in the high temperature phase and 
low temperature phase, respectively, which are taken to 
be ;3i = /3c - 100/iVy, /32 /3c + lOO/TVy, i.e., we ensure 
that the window defining the transition narrows as the 
width of the heat capacity peak narrows. Fig [TO] shows 
the density of the entropy difference AS/Ny as a func- 
tion of Ny, which, rather than remaining constant, is a 
monotonically increasing function. Thus the entropy of 
the system is super-extensive. If the ratio fin/i^L of the 
model scales like Ny as argued above, AS/Ny will 
have the form AS/Ny = 7' \nNy + b. The best fit fine 
using this logarithmic function is also shown in Fig |10[ 
which is consistent with a super-extensive entropy, with 
7' ~ 0.065. 

D. Geometric properties 

In this sub-section, we analyze some geometric proper- 
ties of the two phases: if a geometric property is distinct 
in the two phases, it can serve as an order parameter that 
signals the phase transition. 

As was mentioned before, because the low energy 
graphs are nearly triangulations for our Hamiltonian, 
it is useful to introduce an order parameter that mea- 
sures how similar graphs are to triangulations. For this 
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FIG. 12. Average distance (d) between pairs of vertices, 
plotted as a function of inverse temperature /?, for the sys- 
tem of size Nv = 1000. (d) is first averaged over all pairs 
of vertices in a given snapshot, and then averaged over all 
snapshots at a given temperature. The vertical bars at each 
data point indicate the standard deviation between snapshots: 



purpose we can study the distribution of edge valencies, 
where the edge valency is defined as the number of tri- 
angles that an edge is part of. In a perfect triangulation 
of a surface without boundaries, the edge valencies are 
always two, so we expect that at low temperatures, the 
distribution of edge valency should approximate a delta 
function around two. The distribution of edge valencies 
for the system of size Ny = 1000, X = 100 is shown in 
Figure [TT] as a function of temperature. Indeed, almost 
all edges have edge valency two at temperatures below 
the transition temperature. Near the transition temper- 
ature however, there is a sudden change in the distribu- 
tion of edge valencies: above the transition temperature, 
edge valencies both above and less than two appear. 

Another quantity that is useful as an order parame- 
ter is the average distance between all pairs of vertices, 
denoted by (cZ), where the bar means averaging over 
all pairs of vertices in a graph, and the angle braket 
means averaging over samples of an equilibrium simula- 
tion. We expect that above the phase transition temper- 
ature, graphs will exhibit "small world" topologies and 
thus (d) will be relativity small. The quantity (d) gives 
the characteristic linear size of the graphs. Figure 12 
plots (d) vs. inverse temperature /3, for Ny — 1000. In- 
deed, the low temperature phase has a larger (^J) than the 
high temperature phase; low temperature graphs tend to 
have much more structure than high temperature graphs. 
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FIG. 13. (Color online) The average distance (d) between 
pairs of vertices as a function of the system size A'^v' (discrete 
points), and the best fit lines using a square root function 
(green dashed lines) , using a power function (red solid lines) , 
and using a logarithmic function (black solid lines). Plots are 
shown both above the transition (/? — 1.0) in panel (a), and 
below the transition (/3 — 2.0) in panel (b). For each best fit 
line, its expression and p-value are also shown, where the p- 
values are calculated for the null hypothesis that the residues 
(dfit — {d))/5d come from a normal distribution with vari- 
ance smaller than 1. For both temperatures, the logarithmic 
function gives the best fit to the measured data. 



resulting in larger values of (S^. 

In Figure 13 the average distance (cT) is shown as a 
function of the system size iVy, at /3 = 1.0 (above the 
transition) and at /? = 2.0 (below the transition). The 
best fit lines using a logarithmic function and using a 



power function are also shown in Figure 13 The p- value 
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FIG. 14. (Color online) Log-log plot of (Nr), the thermally 
averaged number of vertices within a distance r, as a function 
of r; the slope gives the dimensionality of the system, which 
in this case is distance-dependent. The plot shown is for the 
system with size iVv = 2000, at /3 = 1.0 (blue) and at /3 = 2.0 
(red) , which bracket the transition. For the same system, the 
inset shows the fractal dimension as a function of r. 

for each best fit line is calculated for the null hypothesis 
that the residues (dfit — {d^)/^d. come from a normal dis- 
tribution with variance smaller than 1, so that a higher 
p-value indicates a better model. These relations be- 
tween ((i) and Ny can be understood by comparing with 
random graphs, which generally display "small- world" 
connectivity, with average distances growing logarithmi- 
cally with the number of vertices [1 . In our model, the 
Hamiltonian only constrains the graphs locally, so these 
graphs satisfy small-world behavior in the high tempera- 
ture phase accurately, as shown by the logarithmic best 
fit line in Figure [iSf^a). For the low temperature phase, 
we can define an effective scaling dimension (see e.g. [25] ) 



D, 



d\nN^ 



V 



d\ll{d)' 



(8) 



On a non-fractal surface, (fi) ^ N^^^, i.e., Ds = 2. 
However, it is seen from Figure [l3jb) that the residu- 
als with the square root function are too large. If we 
take Ds as a parameter in the fitting, a power law func- 
tion with Ds — 3.5 is a much better fit to the empirical 
scaling. Perhaps surprisingly however, the logarithmic 
fimction is still the best fit function, indicating that the 
low temperature graphs still display small-world connec- 
tivity. Enforcing a power-law fit at every system size, i.e. 
{d} ^ Ny^^"^^''\ would induce a variable dimensionality 
in the exponent. 

Another related definition of dimensionality measures 



the increase in number of vertices with distance from a 
given vertex. On a graph, one can pick an arbitrary cen- 
tral vertex, and count how many vertices Nr have dis- 
tance no greater than r from that center. We can then 
average over both all central vertices and over all equi- 
librium configurations at a given temperature, denoting 
the doubly averag ed volume by (Nr). If (Nr) increases 
with r polynomially, the fractal (Haussdorf) dimension 
can be defined as 



Df = 



d\n{Nr) 
din r 



(9) 



In practice the dimension of the graph may itself depend 
on the radius r, so it makes sense to talk rigorously about 
the dimensionality of a graph only if Z?/ is essentially con- 
stant over some range of r. A log-log plot of (iVV) vs. r 



is shown in Figure [14) for Ny = 2000 at /3 = 1.0 and 
/3 = 2.0, where the slope thus gives the dimensionality 
and is shown in the inset. One can see that the effec- 
tive dimension Df is smaller below the transition. Con- 
sistent with the previous analysis using (|8]), there is no 
well-defined dimension for the graphs, which are small- 
world-like. Instead there is an increasing dimensionality 
with increasing length-scale, until boundary effects of the 
system are felt. The dimensionality has values around 2 
for small values of r, because of the local lattice-like struc- 
ture; it is also small for very large values of r, because a 
finite-sized graph must eventually be bounded, at which 
point {Nr') will no longer increase polynomially at large 
r. Table |l] lists the maximal value of Df[r) for systems 
with different sizes, at inverse temperatures P — l.Q and 
(3 = 2.0. As the tables shows, I?/, max increases with Ny, 
which indicates that as Ny increases, there is no univer- 
sal fractal dimensionality that can be approached by the 
graphs. Instead, the graphs are still small-world. 





P = 1.0 


P = 2.0 




= 1000 


3.62 


2.72 


Ny 


= 1500 


3.99 


3.11 


Ny 


= 2000 


4.26 


3.20 



TABLE I. The maximal value of the fractal dimension Df 
as defined in ^ for systems with Ny = 1000, Ny = 1500 
and Ny = 2000, at inverse temperatures /3 = 1.0 above the 
transition and /? = 2.0 below the transition. 

The small-worldness of the low temperature graphs in 
the bulk limit can be viewed as a consequence of the 
graph Hamiltonian in ([5]), which is a sum of local terms. 
The defects in the manifold are also local - in the bulk 
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these have no effect on the large-scale structure of the 
resulting graphs. This is manifested for finite-size graphs 
by the fact that as Ny increases, the topologies of graphs 
become progressively more complicated, see e.g. Fig- 
ure |4f^c) and Figure |4|^d). The manifolds contain numer- 
ous handles and surface intersections, so that a planar 
dimensionality does not adequately describe the system. 
In this sense there is already the signature in the low tem- 
perature phase of the finite system that the bulk system 
is always disordered. 



E. Correlation functions 

Defects in this model such as those shown in Figure [5] 
contain irregularities which make them differ from part 
of a regular lattice. However, regions far away from them 
may not be affected by their existence, i.e., there may be 
no long-range correlation between such defects. In this 
subsection, we define and calculate correlation functions 
between defect pairs. 

Because valency-7 vertices induce defects, we first mea- 
sure the radial correlation function of valency-7 vertices. 
In general, the correlation between two random variables 
X, Y with expected values fxx , fJ'Y and standard devia- 
tions ax , cry is defined as 



corr(X, Y) 



E[{X - t^xW - fly)] 
axcry 



(10) 



where E is the expectation value operator. In our case, 
we take all pairs of vertices with distance d in a graph; 
X is 1 if the first vertex in a pair has valency 7, and 
otherwise, and Y is defined similarly for the second ver- 
tex. Then the correlation function is averaged over all 
equilibrium samples. The result for Ny = 2000, taken at 
inverse temperatures f3 = 1.62 and /3 = 1.65, which are 
marginally below and above /3c respectively, is shown in 
Figure 15 a). When the distance d is very small (d = 1 or 



2), the correlation function deviates from zero, because of 
the local structure of of the defects (see Figure [5]), which 
in this case induces anti-correlation. For intermediate 
values of d (3 < d < 10), the correlation is very small, 
indicating the defects are uncoupled. However, for large 
values of d, the correlation function becomes negative. 
This is because the valency of a vertex, and the distance 
from this vertex to other vertices, are not independent: 
compared with the valency-6 vertices, the valency-7 ver- 
tices tend to have smaller distances to other vertices. For 
example, for Ny = 2000, /3 = 1.62, the mean distance to 
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FIG. 15. (Color online) Radial correlation function defined 



through ( 10 1 of (a) valency-7 vertices and (b) valency-3 edges. 



Correlations are calculated for the system with size Nv = 
2000 at /3 = 1.62, which is in the high temperature phase, 
and at /3 — 1.65, which is the low temperature phase. 



valency 6 vertices is 7.18, while the mean distance to va- 
lency 7 vertices is 7.04. Thus it is less probable to find 
two valency-7 vertices with a large distance, and hence 
they anti-correlate at large distances. The correlation 
function is quite small over a range of d as one might 
anticipate, but the above global effect makes it difficult 
to quantitatively confirm that defects are decoupled at 
large distance. 

As another measure of the correlation between defects, 
we can measure the radial correlation function of valency- 
3 edges, since their existence indicates deviation of the 
graph from a triangulation of surface. For example, ev- 
ery defect in Figure [5] contains valency-3 edges. The dis- 
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tance between two edges is defined by taking the 4 ver- 
tices defining the two edges, and finding the pair of ver- 
tices with the minimum distance between them. Since a 
pair of edges having a common vertex would then have 
a distance of zero, we add one to the above definition of 
edge distance. The results for Ny = 2000, at /3 = 1.62 
and /3 = 1.65, are shown in Figure [iSjb). At small dis- 
tances {d < 3), there exists short range positive cor- 
relation between the valency-3 edges - the mean force 
between them is attractive, due again to the particular 
structure within a given low energy local defect. At large 
distances (d > 14 for /3 = 1.62, d > 17 for /? = 1.65), the 
correlation function becomes negative, because valency- 
3 edges correlate with valency-7 vertices, which in turn 
anti-correlate at large distances for the reasons described 
above. However for a wide range of intermediate dis- 
tances, this correlation function is also nearly zero, indi- 
cating again that the defect attraction is short-ranged. 




(a) 




V. ADDITION OF COULOMB POTENTIAL 

We found above that as the graph size Ny increased to 
infinity, the transition temperature Tc approached zero 
(Figure [s]) . This is apparently a universal property of 
models based on graphs, due to the super-extensive en- 
tropy of the high temperature random phase. Similar 
arguments appear in the theory of phase transitions of 
low dimensional systems |24| . wherein the non-extensive 
energy cost of defect formation is outweighed at any non- 
zero temperature by the (extensive) free energy due to 
translational entropic gain, so long as interactions are 
sufficiently short-ranged. This analogy motivated us to 
introduce a model with long-ranged interactions between 
defects, anticipating that in such a defect-filled system 
incurs super-extensive energetic cost, which may in turn 
result in a non-zero transition temperature. 

Thus, in addition to the original two terms in the 
Hamiltonian ([5]), we introduce a nonlocal Coulomb po- 
tential term to the Hamiltonian, which gives a repulsive 
force between any pair of degree-7 vertices. 



C3 



v,u&V{G),v^u 



d{v, u) 



(11) 



This is one of the simplest non-local Hamiltonian terms 
that one can add to the original Hamiltonian. The 
Coulomb force is chosen to be repulsive, because most 
of the high temperature states are "small world" , in that 
they have smaller average distances than those of low 




(b) 



FIG. 16. (Color online) Sample configurations for the model 
with Coulomb potential in with C3 = 1.0, for the system 
with number of vertices Ny = 200 and number of extra edges 
X = 20, drawn in three dimensions. Panel (a) shows a typical 
configuration in the high temperature phase with j3 — 1.0; 
panel (b) shows a typical configuration in the low temperature 
phase with /3 = 2.0. 



temperature states, so such a Coulomb potential can sup- 
press the appearance of these "small world" graphs. 

We test the effect of addition of this Coulomb term 
by another set of simulations, in which C3 = 1.0. Fig- 
ure [16] shows the sample drawings of graphs with Ny — 
200, A = 20 (a) at high temperature (/3 = 1.0) and (b) at 
low temperature (/? = 2.0). These temperatures bracket 
the heat capacity peak for the system so that the sys- 
tem is in the disordered and ordered phases respectively 



(Fig. 17). Because of the non-locality of H3, simulations 
are much slower in practice than before and smaller sys- 
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FIG. 17. (Color online) Log-log plot of transition tempera- 
tures /3c as a function of system size Ny, for the model with 
local Hamiltonian in Q (drawn as circles, with best fit drawn 
as solid line), and for the model with Coulomb potential in 
(111 with C3 = 1.0 added to the local Hamiltonian (discrete 
points with error bars). The inset shows the rescaled heat 
capacity C/Ny as a function of inverse temperature /3 for 
systems with the Coulomb potential added, and from which 
the values and uncertainties of /?c values are determined. 
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terns are thus employed: simulations are performed for 
Nv = 80,100,120,150,200 and 250, and X = OANy- 



The inset of Figure 17 shows the rescaled heat capacity 
C/Ny as a function of for several system sizes. The 
rescaling factor is now chosen differently than in Fig- 
ure [Tj because the systems with the Coulomb potential 
have maximal heat capacity approximately proportional 
to Ny. From the maximal heat capacity, the inverse 
transition temperature /?c is determined, and is shown in 



Figure 17 (main panel), in comparison with the /3c values 



without the Coulomb potential. 



From the graph drawings in Figure 16 we can see that 
because of the repulsive Coulomb force, both the high 
temperature and low temperature manifold configura- 
tions become rather extended to achieve longer average 
distances between defects. This may also explain why 
the transition temperature does not change very much 
with Ny: The characteristic linear size of the systems 
is much larger when the repulsive Coulomb potential is 
present, which penalizes the increase in complexity that 
was observed for a local Hamiltonian as Ny increased. 
We thus suspect that the entropy would be extensive for 
the long-ranged interaction model. To quantify this, as 



(b) 



FIG. 18. (Color online) Sample configurations for the model 

-1.0, for the sys- 



with Coulomb potential in (111 with C3 



tem of size Ny — 200, X = 20 drawn in three dimensions. 
Panel (a) shows a typical configuration in the high tempera- 
ture phase with /? = 1.0; panel (b) shows a typical configura- 
tion in the low temperature phase with /3 = 2.0. 



a final check we plot the entropy change between disor- 
dered and ordered phases as a function of Ny in the inset 
of Figure 10 where AS* is calculated by Equation ([7|, and 
/3i = 1.0, P2 = 2.0. As opposed to the entropy difference 
in the original model, A5/iVy of this model is approxi- 
mately constant as Ny increases, i.e. the entropy differ- 
ence is no longer super-extensive, rather it is extensive 
or sub-extensive. 



We also simulate the model with an attractive 
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Coulomb potential, in which C3 = —1.0. Figure [Ts] shows 
sample drawings of graphs with Ny = 200, X — 2Q (a) at 
high temperature (/3 = 1.0) and (b) at low temperature 
(/3 = 2.0). The effect of the attractive potential can be 
observed in these samples, in that the valency-7 vertices 
(green and blue dots) are usually located close together. 
In addition, because a local move must involve a valency- 
7 vertex, the configuration cannot evolve in the regions 
composed of purely valency-6 vertices, and thus the sim- 



ulation is inefficient. As can be seen in Figure 18 'b), in 



the region of valency-6 vertices, the configuration does 
not minimize the Hamiltonian (red dots have positive 
contribution to i?2)j ^nd is not a triangulation. Thus 



Figure 18 b) depicts a long-lived meta-stable state on an 
energy landscape of states characteristic of a frustrated 
system [551 US] • Such a model has numerous local min- 
ima with large reconfigurational barriers between them, 
and consequently glassy relaxation dynamics. 



VI. DISCUSSION 

In this paper we have constructed a graph model with 
a local Hamiltonian that simply enforces valency and 
a graph symmetry between the local graph radius and 
diameter. This Hamiltonian gives rise to an emergent 
manifold at low temperature. The one free parameter 
in the model does not appear in the Hamiltonian but 
as an initial condition of the system. This parameter a 
determines the edge to vertex ratio, which is conserved 
for the system and determines the dimensionality of the 
emergent manifold. When a is slightly larger than 6, 
the low temperature solutions have structural properties 
consistent with triangulations of two dimensional sur- 
faces. Because there is no constraint on the global struc- 
ture of the graph however, the low temperature graphs 
can still retain complicated topologies with small-world 
properties, for which the corresponding manifold shows 
handles, self-intersections, and local defects that deviate 
from the manifold in that a higher embedding dimen- 
sion is necessary to represent them. As a general prop- 
erty of the graph model, the high temperature phase has 
an entropy that grows super-extensively with system size 
Nv- This results in a transition temperature of zero in 
the limit Ny — >■ cxd, so that the infinite manifold is al- 
ways disordered at any finite temperature. Aside from 
finite-universe or diverging coupling constraints as pos- 
sible solutions, we implemented long-range interactions 
between vertex defects with repulsive Coulombic poten- 



tial, to energetically penalize the many graph configu- 
rations with defect arrangements consistent with small 
world topologies. In analogy with low-dimensional con- 
densed matter systems, long-range potentials that couple 
defects induce prohibitive energetic cost to configurations 
that would destroy order entropically, and so recover an 
ordered phase at low temperature. Here, we found that 
such potentials result in a nearly constant transition tem- 
peratures as the size of the graph Ny increases. In addi- 
tion, we found that attractive Coulombic potentials re- 
sult in long-lived metastable states in the simulations. 

Another interesting feature of the model is that the 
low lying energy levels, including the ground state level, 
have large configurational degeneracy. This residual en- 
tropy is due to local defects that can "absorb" extra edges 
without energetic cost. As well, the simulation dynamics 
indicates that the energy barriers between different low 
energy states are not high. Thus at temperatures below 
the phase transition, the degrees of freedom in the sys- 
tem arising from this residual entropy are not frozen. The 
small or zero-energy barriers between degenerate states 
make the low temperature graph system similar to the 
spin ices observed in the spinel structures and pyrochlore 
lattices 



It is intriguing to interpret the low temperature config- 
urations of this graph model as an emergent spacetime - 
a notion other researchers have explored for similar graph 
models [3l IH IMl4l l30] . In this picture, general relativ- 
ity is an effective "hydrodynamic" theory emerging from 
the collective dynamics of more fundamental degrees of 
freedom. The graph model is appealing in that both 
spacetime manifolds and locality emerge in the low tem- 
perature regime of a discrete structure. The graph model 
introduced here gives rise to real, positive distances, so 
the emergent manifold can only be a Euclidean version 
of spacetime. However, the fact that the topology of the 
corresponding surfaces is not preserved even at low tem- 
perature poses a problem in mapping such a model to a 
classical theory of gravity. 

It is still possible that the set of all possible low energy 
graphs in this model could be identified with the phase 
space of a gravity theory before imposing the equation of 
motion, i.e., the space of all possible metrics modulo dif- 
feomorphisms. It would be very interesting to investigate 
whether one could add terms to the action corresponding 
to our graph model, for example those in dynamical trian- 
gulation theory, which reduce to the gravitational action 
in the continuous limit. Such a construction could give a 
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model of emergent gravity from graphs. If this procedure 
were successful, we would expect the fractal dimension of 
the graphs to reduce, simplifying the topology and corre- 
sponding manifold, leading to a smaller contribution to 
the gravitational action. Discrete models of gravity are 
known to have large entropy for the configurations with 
fractal structures however (see, e.g., [311 [32]), so further 
work is needed to understand the energetic conditions 
that would result in a realistic emergent spacetime. 

The graph model in the present study possesses a first 
order phase transition. It is possible that a combined 
model such as those above can have a second-order phase 
transition, or have two successive phase transitions, in 
which one is the current first-order phase transition, and 
the other new phase transition is second-order, corre- 
sponding to the transition between complicated, fractal 
surface-like graphs and extended surface-like graphs. 

The order of the transition, and its potential relevance 
to universality or independence of underlying lattice 
specifics, is a non-issue for the investigation of ordered 
phases below the transition, where correlation lengths 
are finite. Power law correlations calculated in causal 
dynamical triangulation are between graphical elements 
analogous to graviton fields, so that graviton coupling is 
power law as in the classical limit. Space to time ratios 
of simplices have second order transition in this model, 
while the transition involving gravitational coupling is 
first order [33]. In any event, a lattice- like system at a 
critical point would have wildly fluctuating connectivity 
and resemble more a fractal mix of ordered and disor- 
dered states, which is not consistent with an emergent 
manifold. The issue of the universality classes and cor- 
responding exponents of a transition is a separate one 
from the properties of an emergent manifold as a low 
temperature phase below a phase transition. 

In this context however, it is intriguing to speculate on 
the utility of a such a graph theoretical transition to de- 
scribe a transition involving non-local to local causality, 
as might occur in a pre-inflationary universe. Such mod- 
els may address the low-entropy initial condition prob- 
lems that occur in inflationary models [341 135] 
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Appendix A: The weighted histogram analysis 
method 



The weighted histogram analysis method (WHAM) is a 
method to combine the samples from several Monte Carlo 
simulations taken under conditions of different tempera- 
ture and added potential. We employ WHAM to gen- 
erate optimal estimates of energy distributions of the 
graph model. In this model, the energy takes only in- 
teger values between and M = 1.5Nv- Assume that 
S simulations are performed (in our cases, S — A for 
Nv = 1000, 5 = 10 for Ny = 2000), with inverse tem- 
perature Pi, and biasing potential Vi{E), i.e., in the i-th 
simulation, the system is sampled with energy distribu- 
tion n{E) eyi-p{- I3,{E + V^{E))), where n{E) is the yet- 
unknown number of states with energy E. The inverse 
temperature /J^'s are taken to be near the inverse tran- 
sition temperature. Because there is a large free energy 
barrier between the low and high energy phases near the 
transition temperature, a biasing potential is used to ob- 
tain better sampling in the barrier region. The form of 
the biasing potential is taken to be parabolic: 



V^iE) = { 



h , 



0, 
I 0, 



E\<E< E'^, 



E < El 
E > Et 



where the parameters u;, E\ and E^ are chosen by trial- 
and-error to make the energy distribution of each simu- 
lation as flat as possible. 

After performing the simulations, let ni{E) be the 
number of counts of energy E from the i-th simulation, 
and Ni the total number of samples from the i-th sim- 
ulation. From this information, the optimal estimate of 
the probability p^{E) of energy level E at inverse tem- 
perature /3° without any biasing potential is given by 



p\E) = 



J:lln^{E) 
ELiVJ,c,(i?)^ 



(Al) 



where Ci{E) is the biasing factor Ci{E) = exp[—{(3i — 
/3^)E — [3iV{E)], and fi is a normalization constant sat- 
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FIG. 19. For the system of size Nv ~ 1000, the acceptance 
ratio of Monte Carlo moves in the simulations as a function 
of inverse temperature /3. 



isfying 



M 



/-i = 5]c.(i?)/(£;). 



(A2) 



To solve these equations, we take an arbitrary set of ini- 
tial values for fi (namely /f = 1), and apply (Al) and 



(A2) iteratively to find the solution to these equations. 
After finding it is then straightforward to calculate 

the average energy and heat capacity at inverse temper- 
ature 



Appendix B: Acceptance ratio 

As a practical matter, we plot the acceptance ratio as 
a function of (5 for Ny = 1000 in Figure [19] The low 
energy phase occupied at large values of f3 has a much 
lower acceptance ratio than the high energy phase, both 
because of the lower temperature and because the low en- 
ergy graphs have much more structural constraints, and 
thus have more rigidity with respect to the local moves. 
However, because some the local defects cost little or no 
energy, low energy graphs still have non-zero acceptance 
ratio, and so are able to undergo dynamics during the 
simulations. 
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